I'm gonna overwrite a lot of this notebook's old content. I changed the way I'm calculating wt, and wanna test that my training worked.



In [53]:

    
from pearce.emulator import *
from pearce.mocks import cat_dict
import numpy as np
from os import path



In [54]:

    
import matplotlib
#matplotlib.use('Agg')
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()



In [55]:

    
from GPy.models import GPKroneckerGaussianRegression, GPRegression
from GPy.kern import *
import h5py



In [56]:

    
training_file = '/home/users/swmclau2/scratch/xi_gg_zheng07_cosmo_v3/PearceXiggCosmo.hdf5'
test_file = '/home/users/swmclau2/scratch/xi_gg_zheng07_cosmo_test_v3/PearceXiggCosmoTest.hdf5'



In [57]:

    
f = h5py.File(training_file, 'r')



In [58]:

    
Ys = []
rbin = 0
for i in xrange(40):
    Ys.append(f['cosmo_no_%02d'%i]['a_1.000']['obs'].value[:, rbin])



In [59]:

    
n_hods = 100
start_idx = 0
X1 = f.attrs['cosmo_param_vals']
X2 = f.attrs['hod_param_vals'][start_idx:start_idx+n_hods]
Y = np.vstack(Ys)[:, start_idx:start_idx+n_hods]



In [60]:

    
X1[-10:,:]









    Out[60]:





array([[  2.19279200e-02,   1.23913700e-01,  -1.22368100e+00,
          9.55235800e-01,   3.11828200e+00,   7.34315200e+01,
          4.06625000e+00],
       [  2.12396400e-02,   1.27578700e-01,  -1.38164500e+00,
          9.56072500e-01,   3.07641300e+00,   7.37599700e+01,
          3.34375000e+00],
       [  2.24968800e-02,   1.12845000e-01,  -9.25757800e-01,
          9.49484200e-01,   3.04281000e+00,   6.84034200e+01,
          3.98125000e+00],
       [  2.33731900e-02,   1.14970300e-01,  -8.75220700e-01,
          9.89173100e-01,   3.14916300e+00,   6.60523100e+01,
          3.64125000e+00],
       [  2.27557300e-02,   1.22229300e-01,  -1.03229000e+00,
          9.50028100e-01,   3.10694000e+00,   6.90703800e+01,
          3.13125000e+00],
       [  2.34107700e-02,   1.07640600e-01,  -6.12819700e-01,
          9.95567600e-01,   3.14028900e+00,   6.16947200e+01,
          3.04625000e+00],
       [  2.19989400e-02,   1.21296700e-01,  -1.10822200e+00,
          9.67361000e-01,   3.17942400e+00,   7.04112900e+01,
          3.30125000e+00],
       [  2.28744700e-02,   1.09744000e-01,  -8.48745500e-01,
          9.77622700e-01,   3.07228700e+00,   6.67262000e+01,
          3.51375000e+00],
       [  2.37123900e-02,   1.14974700e-01,  -9.55007500e-01,
          9.76590700e-01,   3.05350800e+00,   6.97468000e+01,
          4.10875000e+00],
       [  2.17407200e-02,   1.20070000e-01,  -9.41052900e-01,
          9.60207200e-01,   3.09270400e+00,   6.47043000e+01,
          4.19375000e+00]])



In [61]:

    
X2[-10:,:]









    Out[61]:





array([[ 12.96736737,   0.06126126,  14.3966967 ,   1.26756757],
       [ 13.04444444,   0.4027027 ,  14.18548549,   1.10960961],
       [ 13.52692693,   0.18558559,  14.5028028 ,   1.22312312],
       [ 12.82122122,   0.27072072,  13.73403403,   1.11021021],
       [ 12.91331331,   0.4018018 ,  14.44374374,   0.7960961 ],
       [ 12.97337337,   0.2518018 ,  14.55085085,   1.1006006 ],
       [ 13.52892893,   0.29279279,  13.85915916,   1.04054054],
       [ 12.83423423,   0.05720721,  13.89119119,   0.98108108],
       [ 13.41981982,   0.30675676,  14.11341341,   1.25435435],
       [ 13.58998999,   0.32162162,  13.8981982 ,   0.74624625]])



In [62]:

    
X1.shape, X2.shape, Y.shape









    Out[62]:





((40, 7), (100, 4), (40, 100))



In [63]:

    
# how to add training errors?



In [64]:

    
c = np.eye(X1.shape[0])
print c.shape
print c.diagonal()









    



(40, 40)
[ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  1.  1.  1.  1.]

K1 = RBF(input_dim = 7, ARD = False)+ Fixed(input_dim = 7, covariance_matrix=np.eye(X1.shape[0]))#+Bias(input_dim=7) + Linear(input_dim = 7, ARD = True)# + Bias(input_dim=7)# + White(input_dim=7) K2 = RBF(input_dim=4, ARD = False)#+ Linear(input_dim = 4, ARD = False) + Bias(input_dim=4)# + White(input_dim=4)



In [65]:

    
K1 = RBF(input_dim = 7, ARD=False)#+ Linear(input_dim = 7, ARD = False) + Bias(input_dim=7)# + White(input_dim=7)
K2 =  RBF(input_dim=4, ARD = False)#Linear(input_dim = 4, ARD = True) + Bias(input_dim=4)# + White(input_dim=4)



In [66]:

    
model = GPKroneckerGaussianRegression(X1, X2, Y,K1, K2 )



In [67]:

    
model.optimize_restarts(num_restarts=10)









    



Optimization restart 1/10, f = -13151.9502508
Optimization restart 2/10, f = -13151.9502413
Optimization restart 3/10, f = -13151.9502523
Optimization restart 4/10, f = -13151.9502642
Optimization restart 5/10, f = -13151.9502335
Optimization restart 6/10, f = -13151.950267
Optimization restart 7/10, f = -13151.9502604
Optimization restart 8/10, f = -13151.9502683
Optimization restart 9/10, f = -13151.9502764
Optimization restart 10/10, f = -13151.9502728






    Out[67]:





[<paramz.optimization.optimization.opt_lbfgsb at 0x7fa40165ee50>,
 <paramz.optimization.optimization.opt_lbfgsb at 0x7fa401615bd0>,
 <paramz.optimization.optimization.opt_lbfgsb at 0x7fa4015ed290>,
 <paramz.optimization.optimization.opt_lbfgsb at 0x7fa4015ed210>,
 <paramz.optimization.optimization.opt_lbfgsb at 0x7fa401651e10>,
 <paramz.optimization.optimization.opt_lbfgsb at 0x7fa401651f50>,
 <paramz.optimization.optimization.opt_lbfgsb at 0x7fa401651fd0>,
 <paramz.optimization.optimization.opt_lbfgsb at 0x7fa401651ed0>,
 <paramz.optimization.optimization.opt_lbfgsb at 0x7fa401651f10>,
 <paramz.optimization.optimization.opt_lbfgsb at 0x7fa401651dd0>]



In [85]:

    
model.kern1









    Out[85]:





  rbf.               value constraints priors
  variance   69.5743322297     +ve          
  lengthscale 6.14360809854     +ve



In [68]:

    
K1.param_array









    Out[68]:





array([ 69.57433223,   6.1436081 ])



In [69]:

    
K2.param_array









    Out[69]:





array([ 45.88148894,   1.97012146])



In [70]:

    
f2 = h5py.File(test_file, 'r')



In [71]:

    
Y2s = []
for i in xrange(35):
    Y2s.append(f2['cosmo_no_%02d'%i]['a_1.000']['obs'].value[:, rbin])



In [72]:

    
testX1 = f2.attrs['cosmo_param_vals']
testX2 = f2.attrs['hod_param_vals']
testY = np.vstack(Y2s)



In [73]:

    
print testX1.shape
print testX2.shape
print testY.shape









    



(35, 7)
(1000, 4)
(35, 1000)



In [74]:

    
print testX1[0,:]
print testX2[0,:]









    



[  2.32629000e-02   1.07830000e-01  -7.26513000e-01   9.80515000e-01
   3.03895000e+00   6.32317000e+01   2.95000000e+00]
[ 13.48388388   0.26666667  14.14344344   1.26636637]



In [75]:

    
predY, _ = model.predict(testX1, testX2)



In [76]:

    
print predY.shape



In [77]:

    
print predY[0]









    



[ 3.81625036]



In [78]:

    
np.median( np.abs( (10**predY[:,0] - 10**testY.flatten(order='F'))/10**testY.flatten(order='F') )  )









    Out[78]:





0.029989317626029389



In [79]:

    
np.mean( np.abs( (10**predY[:,0] - 10**testY.flatten(order='F'))/10**testY.flatten(order='F') )  )









    Out[79]:





0.037938955447990814



In [80]:

    
predY[:,0]









    Out[80]:





array([ 3.81625036,  3.81625036,  3.81625036, ...,  3.8831422 ,
        3.8831422 ,  3.8831422 ])



In [81]:

    
testY.flatten(order='F')









    Out[81]:





array([ 3.75435156,  3.76497942,  3.76070541, ...,  3.89549228,
        3.89280012,  3.89538578])



In [82]:

    
k_dict = K1.to_dict()



In [83]:

    
k_dict['parts'].keys()









    




KeyErrorTraceback (most recent call last)
<ipython-input-83-4a2605541446> in <module>()
----> 1 k_dict['parts'].keys()

KeyError: 'parts'



In [ ]: